A deep dive into the various mathematical notation systems and R packages for fitting multilevel hierarchical models
This is a blog post that summarizes the things I found helpful while going through the various mathematical notation systems and R syntax for multilevel models. Some quick takeaways:
R (be it lm, lmer, glm, glmer, or brms), the R formula most closely resembles the composite notationequitiomatic, the quickest way to understand the math (for all fitted models, frequentist and bayesian) is to probably see the equation for yourself using the extract_eq function.1 See Kenneth Iverson’s writing for more details.
For this post we will consider main effects only (i.e., no interaction terms).
McElreath (2020) page 415:
Multilevel models go by many different names, and some statisticians use the same names for different specialized variants, while other use them all interchangeably. The most common synonymns for “multilevel” are HIERARCHICAL and MIXED EFFECTS. The type of parameters that appear in multilevel models are most commonly known as RANDOM EFFECTS, which itself can mean very different things to different analysts and in different contexts. And even the innocent term “level” can mean differeng things to different people. There’s really no cure for this swamp of vocabulary aside from demanding a mathematical or algorithmic definition of the model. Otherwise, there will always be ambiguity.
McElreath (2020) footnote 188:
I adopt the terminology of Gelman (2005), who argues that the common term random effects hardly aids with understanding, for most people. Indeed, it seems to encourage misunderstanding, partly becasue the terms fixed and random mean different things to different statisticians. See pages 20 – 21 of Gelman’s paper. I fully realize, however, that by trying to spread Gelman’s alternative jargon, I am essentially spitting into a very strong wind.
Often times we encounter data with some hierarchical nature, e.g., students’ scores from schools from different districts, etc. Here we will use the data set used by this chapter by Stef van Buuren in Flexible Imputation of Missing Data. This is data on school performance of grade 8 pupils in Dutch schools. There data were collected by Brandsma and Knuver (1989) and were used as the primary examples in Chapters 4 and 5 of Snijders and Bosker (2011). The data contains a mix of both individual/pupil-level measurements and group/school-level measurements. We will concentrate on four variables:
sch: school numberlpo: language test post, outcome at pupil level.sex: sex of pupil (coded 0 - 1)den: school denomination, predictor at school level (coded 1 - 4)data("brandsma", package = "mice")
df <- brandsma %>% select(sch, lpo, sex, den)
# there is NA data so we filter those out
df <- df %>% na.omit() %>%
# turn the other columns into factors
mutate(sch = as.factor(sch),
sex = as.factor(sex),
den = as.factor(den))
# df %>% glimpse()
# df %>% summary()
df %>% head(5) sch lpo sex den
2 1 50 1 1
3 1 46 0 1
4 1 45 0 1
5 1 33 0 1
6 1 46 0 1df %>% ggplot(aes(x = lpo, color = sex)) +
geom_density()
df %>% ggplot(aes(x = lpo, color = den)) +
geom_density() +
facet_wrap(~sex)
We want to create a model that predicts whether the students’ performance (lpo) is influenced by their (sex) and which school they’re from. We will start with …
3 or potentailly more
\[\begin{align} \texttt{lpo}_{ij} &= \beta_{0j} + \beta_{1j}\texttt{sex}_{ij} + \epsilon_{ij} \\ \beta_{0j} &= \gamma_{00} + \gamma_{01} \texttt{den}_j + u_{0j} \\ \beta_{1j} &= \gamma_{10} \end{align} \]
where
nlme, lmer)By substituting the level-2 equations into the level-1 equation, we get the following one-line equation:
\[ \texttt{lpo}_{ij} = \beta_0 + \beta_1 \texttt{sex}_{ij} + \beta_2 \texttt{den}_j + u_{0j} + \epsilon_{ij} \]
It’s usually assumed that \(u_{0j} \sim \mathcal{N}(0, \sigma_u^2)\) and \(\epsilon_{ij} \sim \mathcal{N}(0, \sigma_\epsilon^2)\). Note that \(\beta_0 = \gamma_{00} + \gamma_{10}\) from the previous model.
m2 <- lmer(lpo ~ 1 + sex + den + (1 | sch), data = df)
summary(m2)Linear mixed model fit by REML ['lmerMod']
Formula: lpo ~ 1 + sex + den + (1 | sch)
Data: df
REML criterion at convergence: 25717.2
Scaled residuals:
Min 1Q Median 3Q Max
-4.0026 -0.6520 0.0891 0.7292 2.6420
Random effects:
Groups Name Variance Std.Dev.
sch (Intercept) 15.70 3.962
Residual 61.75 7.858
Number of obs: 3648, groups: sch, 198
Fixed effects:
Estimate Std. Error t value
(Intercept) 38.0682 0.5843 65.157
sex1 2.2725 0.2662 8.536
den2 3.2756 0.7697 4.255
den3 1.1432 0.8273 1.382
den4 4.5931 1.5181 3.026
Correlation of Fixed Effects:
(Intr) sex1 den2 den3
sex1 -0.220
den2 -0.722 -0.001
den3 -0.673 0.005 0.510
den4 -0.368 0.007 0.278 0.259For a quick double check:
extract_eq(m2)\[ \begin{aligned} \operatorname{lpo}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(\operatorname{sex}_{\operatorname{1}}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\gamma_{0}^{\alpha} + \gamma_{1}^{\alpha}(\operatorname{den}_{\operatorname{2}}) + \gamma_{2}^{\alpha}(\operatorname{den}_{\operatorname{3}}) + \gamma_{3}^{\alpha}(\operatorname{den}_{\operatorname{4}}), \sigma^2_{\alpha_{j}} \right) \text{, for sch j = 1,} \dots \text{,J} \end{aligned} \]
# version that uses coefficient
extract_eq(m2, use_coefs = TRUE)\[ \begin{aligned} \operatorname{\widehat{lpo}}_{i} &\sim N \left(38.07_{\alpha_{j[i]}} + 2.27_{\beta_{1}}(\operatorname{sex}_{\operatorname{1}}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(3.28_{\gamma_{1}^{\alpha}}(\operatorname{den}_{\operatorname{2}}) + 1.14_{\gamma_{2}^{\alpha}}(\operatorname{den}_{\operatorname{3}}) + 4.59_{\gamma_{3}^{\alpha}}(\operatorname{den}_{\operatorname{4}}), 3.96 \right) \text{, for sch j = 1,} \dots \text{,J} \end{aligned} \]
Note that the above notation is the Gelman and Hill (2007) approach, where \(j[i]\) indicates group membership of individual \(i\).
The most common notation is as follows: \[\mathbf{y} = X \mathbf{\beta} + Z \mathbf{u} + \mathbf{\epsilon}\] where
Still the same model can be written in matrix notation, just as we have seen before with the linear model. A linear mixed model (LMM) is typically written as:
\[ y_j = X_j \beta + Z_j u_j + \epsilon_j \]
where \(y_j = \texttt{lpo}_j\) is a column vector containing the scores in cluster \(j\), where \(X_j = (1, \texttt{den}_j, \texttt{sec}_j)\) is the design matrix.
There are some nice introductions on this blog.
4 ”these components” are ”… variance components, which are the parameters governing the stochastic terms in a model”
In this blog post by Solomon Kurz you get a mixture of components, but he mostly specify things in the style of McElreath.
brms
If the model is not necessarily linear or Gaussian, we could use glm instead of lm and glmer instead of lmer.
Let’s say you don’t know what “unconditional means model” mean5, but you see its mathematical equation as the following: \[y_{ij} = \beta_0 + u_{0j} + \epsilon_{ij}\]
5 Here’s more reading on unconditional mean model
The model posits that the score of a randomly selected student \(i\) in a randomly selected school \(j\) can be represented as a sum of:
(I don’t think the above language is the most correct description … it seems very forced … )
m3 <- lmer(lpo ~ 1 + (1| sch), data = df)
m3 Linear mixed model fit by REML ['lmerMod']
Formula: lpo ~ 1 + (1 | sch)
Data: df
REML criterion at convergence: 25815.55
Random effects:
Groups Name Std.Dev.
sch (Intercept) 4.239
Residual 7.936
Number of obs: 3648, groups: sch, 198
Fixed Effects:
(Intercept)
40.9 extract_eq(m3)\[ \begin{aligned} \operatorname{lpo}_{i} &\sim N \left(\alpha_{j[i]}, \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for sch j = 1,} \dots \text{,J} \end{aligned} \]
“Multilevel models contain two forms of stochastic variation, usually called variance components, captured in random effects (\(u_{0j}\)) and error (\(\epsilon_{ij}\) ) terms, the former representing systematic differences between groups and the latter unexplained, subject-level differences within groups.”
Let’s suppose you’re a psycholinguist. You administed some tests to a group of participants, where each participant answerer different item types multiple types. The model is not necessarily “nested” but given that a same participant is sampled multiple times, you can group by individual participants (Subject) and you can also group by item types (Item). This is the case demonstrated in Barr et al. (2013). Note that \(X_i\) is just an indicator variable, of whether the item is type A versus type B (? to be confirmed). Barr et al. (2013) provided a handy table for models and their associated lmer syntax, which I will recreate:
| math | lm/lmer formula |
|---|---|
| \(Y_{si} = \beta_0 + \beta_1 X_i + e_{si}\) | lm(Y ~ 1 + X) |
| \(Y_{si} = \beta_0 + S_{0s} + \beta_1 X_i + e_{si}\) | Y ~ 1 + X + (1 | Subject) |
| \(Y_{si} = \beta_0 + S_{0s} + (\beta_1 + S_{1s}) X_i + e_{si}\) | Y ~ 1 + X + (1 + X | Subject) |
| \(Y_{si} = \beta_0 + S_{0s} + I_{0i} + (\beta_1 + S_{1s}) X_i + e_{si}\) | Y ~ 1 + X + (1 + X | Subject) + (1 | Item) |
| \(Y_{si} = \beta_0 + S_{0s} + \beta_1 X_i + e_{si}\) | Y ~ 1 + X + (1 | Subject) + (1 | Item) |
| \(Y_{si} = \beta_0 + I_{0i} + (\beta_1 + S_{1s}) X_i + e_{si}\) | Y ~ 1 + X + (0 + X | Subject) + (1 | Item) |
Note that the above math notation is not the most conventional. Everything except the last part I don’t quite get.
My mental model of how hierarchical models work is heavily influenced by McElreath (2020). In his book, McElreath said:
I adopt the terminology of Gelman (2005), who argues that the common term random effects hardly aids with understanding, for most people. Indeed, it seems to encourage misunderstanding, partly because the terms fixed and random effects mean different things to different statisticians. See pages 20–21 of Gelman’s paper. I fully realize, however, that by trying to spread Gelman’s alternative jargon, I am essentially splitting into a very strong wind.
The TLDR is that there have been a lot of different definitions of fixed and random effects in different contexts, across different disciplines. Gelman chose to sidestep the overloaded terms and define
… effects (or coefficients) in a multilevel model as constant if they are identical for all groups in a population and varying if they are allowed to differ from group to group.
What this implies in terms of how the math notation6 is that
6 Referring to section 12.1 Notation from Gelman and Hill (2007)
According to McElreath (2020), it is a compact form of notation, that are intimately related to a “common matrix algebra description of linear models”.
Xiaoying Pu provided a cheat sheet for reading/writing R formulas used in lmer and brms.
In Barr et al. (2013), the provide a nice table.
In Ga& et al. (2013), they looked at the software for fitting multilevel models, and they had a nice table7 about statistical software packages containing modules/functions designed to fix various classes of mixed-effects models. Below I will reproduce this table, but focusing only on the R side of things:
7 Table 26.1
R package |
Linear Mixed Models LMM |
Generalized Linear Mixed Models GLMM |
Nonlinear Mixed Models NLMM |
|---|---|---|---|
nlme |
lme (1999) |
nlme (1999) |
|
lme4 |
lmer (2005) |
glmer (2006) |
nlmer (2007) |
brms |
brm (2017) |
brm and specify family = bernoulli() or some other non-Gaussian family |
brm and specify nl = TRUE |
See also the footnote of the wikipedia article for multi-level models.